Lanczos algorithm with Matrix Product States for dynamical correlation functions 
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The density-matrix renormalization group (DMRG) algorithm can be adapted to the calculation 
of dynamical correlation functions in various ways which all represent compromises between com- 
putational efficiency and physical accuracy. In this paper we reconsider the oldest approach based 
on a suitable Lanczos-generated approximate basis and implement it using matrix product states 
(MPS) for the representation of the basis states. The direct use of matrix product states combined 
with an ex-post reorthogonalization method allows to avoid several shortcomings of the original ap- 
proach, namely the multi-targeting and the approximate representation of the Hamiltonian inherent 
in earlier Lanczos-method implementations in the DMRG framework, and to deal with the ghost 
problem of Lanczos methods, leading to a much better convergence of the spectral weights and 
poles. We present results for the dynamic spin structure factor of the spin-1/2 antiferromagnetic 
Heisenberg chain. A comparison to Bethe ansatz results in the thermodynamic limit reveals that the 
MPS-based Lanczos approach is much more accurate than earlier approaches at minor additional 
numerical cost. 

PACS numbers: 71.10.Pm,71.10.Fd, 78.20 Bh 



I. INTRODUCTION 

Dynamical correlation functions. A central question 
in physics is the calculation of time-dependent correla- 
tion functions of the type iGAit',t) = {0\A^t')A{t)\0), 
where A is some operator, |0) the ground state and t' > t 
for causality. In the case of a system given by a time- 
independent Hamiltonian iJ, it is often more convenient 
to Fourier transform to frequency space, arriving at some 
dynamical correlation function Ga(j^ + i?y), taking the 
form 



GA(w + i77) = (0|# 



1 



eo + uj + ir] - H 



A\0), (1) 



where eq is the ground state energy and rj an infinites- 
imal positive number introduced for convergence. The 
precise calculation of such and similar dynamical corre- 
lation functions is highly important in condensed matter 
physics as they can be directly compared to experiments 
measuring the response of a many-body system to an ex- 
ternal perturbation, for example in angle-resolved pho- 
toemission spectroscopy or inelastic neutron scattering, 
upon suitable choices of A. From a methodological point 
of view, dynamical correlation functions are important 
ingredients for several approximation schemes in con- 
densed matter theory, for example cluster perturbation 
theory (CPT), variational cluster approximation (VGA), 
dynamical mean field theory (DMFT) and others (see 
Ref. [l| for a comprehensive review). In these approaches, 
one needs so-called impurity or cluster solvers, which pro- 
vide results for Green's functions entering the further 
many-body calculations of the method. Usually, analyti- 
cal solutions arc not available and one has to resort to nu- 



merical approaches; dynamical correlation functions have 
been calculated in numerical methods as diverse as the 
numerical renormalization group, exact diagonalization, 
quantum Monte Garlo and the density-matrix renormal- 
ization group (DMRG). The latter is at the basis of this 
work. 

Dynamical correlation functions in DMRG. The 
density-matrix renormalization group (DMRG)^"— has 
established itself as the most powerful method for the 
calculation of ground states of strongly correlated sys- 
tems in one dimension, such that its extension to dynam- 
ical quantities such as Ga(w -f- ir?) has been an obvious 
focus of research. However, the calculation of dynami- 
cal quantities within the DMRG is still a difficult task, 
for which various approaches have been developed which 
show different combinations of advantages and disadvan- 
tages, which we will briefly outline. 

Hallberg^ made a first attempt to tackle this task by 
using a continued fraction expansion of dynamical cor- 
relation functions based on a Lanczos algorithm which 
had already been successfully employed in the context 
of exact diagonalization calculations j^i^ This algorithm 
is fast, easily implemented in a ground state DMRG 
program, and resolves correlation functions Ga (w -I- ir]) 
whose spectral weight is concentrated in relatively sharp 
cj-peaks very well. Its drawback is that it does not re- 
solve correlators with broadly distributed spectral weight 
very well. However, it has found numerous successful 
applications 1^"— Recently, some of the authors have pro- 
posed an adaptive Lanczos methodic which improves 
over the original method of Hallberg. 

Shortly after Hallberg's original proposal, Ramasesha 
et al^ and later Kiihncr and Whitcl- introduced the cor- 



rection vector method, which since then has successfully 
been applied to many model systems.—^— The correction 
vector method achieves high precision, albeit at much 
higher numerical cost than the Lanczos method. Fourier- 
transforming time-dependent DMRG data offers another 
approach to dynamical spectral functions.—^— However, 
accessing long time scales, to obtain good frequency res- 
olution, is limited by a rapid increase of entanglement 
and the resulting increasing demand of computational re- 
sources. Very recently, dynamical spectral functions were 
also calculated using an expansion of spectral functions 
in terms of Chebychev polynomials^^ While early results 
look very promising, the full potential of this method is 
still unexplored. 

DMRG in the MPS framework. In an independent de- 
velopment, it has been recognized early in the history of 
DMRG that DMRG produces quantum states that take 
the very specific form of so-called matrix product states 
(MPS).— Indeed, DMRG - with minor variations to the 
original formulation - can be recognized to be simply a 
variational search for the lowest energy state within that 
state class i^^i^^ In the context of ground state DMRG this 
insight allows a more elegant formulation of the method 
and highlights connections to entanglement theory, but 
does not yield a more performing algorithmi^ This story 
changes, however, once one is confronted with an algo- 
rithmic situation where more than one quantum state of 
interest is being calculated, not just say the ground state. 
In the original DMRG, the problem of an acceptable ap- 
proximation for several states at one time is solved by the 
procedure of multi-targeting: the reduced density oper- 
ator that governs the choice of reduced bases in DMRG 
is set up from a weighted sum of reduced density opera- 
tors formed for each state individually. This implies that 
none of the states will be represented with the accuracy 
the method could have achieved for it alone. In the MPS 
formulation, each state carries its own implicit (optimal) 
choice of reduced basis, hence is approximated more pre- 
cisely; moreover, the formalism automatically deals with 
the fact that each state has been approximated differ- 
ently. Another advantage is that upon the introduction of 
matrix product operators (MPOs) the Hamiltonian can 
be expressed exactly whereas in DMRG it is expressed 
in the reduced basis chosen by DMRG to represent the 
state(s) of interest. 

Combining Lanczos and MPS. Given the range of op- 
tions for calculating dynamical correlation functions in 
the DMRG framework, it would be highly desirable if the 
Lanczos method, which is computationally very unde- 
manding and has very few fine-tuning parameters, could 
be improved in its range of applicability while maintain- 
ing low computational cost. The purpose of this paper 
is to provide a substantial improvement of the current 
Lanczos-based methods. This will be achieved by mak- 
ing use of the methodological progress due to the use of 
MPS and MPOs and a few additional modifications: 

(i) As already noted in Ref. [ij, the MPS formulation 
should be much more suited than original DMRG for the 



Lanczos method, because it generates a long sequence 
of quantum states that call for optimal representation 
beyond multi-targeting; moreover, they are obtained by 
subsequent applications of the Hamiltonian, such that 
an exact representation of the Hamiltonian should be 
very useful. In this work we will therefore use the MPS 
formulation^!^"— of DMRG to implement the Lanczos 
algorithm. Every Lanczos state will be represented by 
one MPS and the only approximation comes from the 
compression of the state. 

(ii) The problem of numerical loss of orthogonality be- 
tween Lanczos states is at first sight exacerbated in the 
MPS formulation due to the required, but inexact com- 
pression of MPS after addition of MPSs or application of 
the Hamiltonian operator to an MPS (both inevitable in 
the Lanczos algorithm) in comparison to the exact Lanc- 
zos method. It can, however, be removed by an ex-post 
reorthogonalization method. 

(iii) The method allows direct access to spectral poles 
and weights without broadening, and we can do finite- 
size extrapolations to the thermodynamic limit using a 
special rescaling procedure. 

The new approach combines the advantages of a very 
good representation of individual states, largely eliminat- 
ing the so-called "ghost" problem of Lanczos approaches, 
and a removal of the broadening issue. 

Outline of the paper. In Section |TT1 we will introduce 
the model to be used for the exposition of the new method 
and as a test case, the spin-i Heisenberg model. In Sec- 
tion [Hll we will review the basic Lanczos algorithm for 
dynamical correlation functions, and its inherent prob- 
lems (finite-size broadening, loss of orthogonality) and 
briefly discuss the status of the current Lanczos-based 
methods, before we introduce the new approach (Section 
lIVp . In order to test the performance of the algorithm 
we have calculated the dynamical structure factor of the 
antiferromagnetic spin-i Heisenberg chain (Section |V|. 
This quantity is not characterized by a single magnon 
band, where the pole structure is very simple, but rather 
a multi-spinon continuum. Such complex structures are 
generally hard to access by the Lanczos method, and 
hence provide a relevant test case to see how this ap- 
proach performs for systems with non-trivial dynamical 
properties. To this end we compare our results with those 
from the thermodynamic Bethe ansatz. Wc find that the 
MPS-based approach leads to much more accurate results 
than earlier Lanczos methods. The paper is concluded by 
a discussion (Section IVI[1 . 



II. MODEL 



As model to test our algorithm we choose the antifer- 
romagnetic spin-i Heisenberg chain with Hamiltonian 



ff=^S,S 



i+l 



Here Si denote the usual spin-operators at site i and we 
have set the antiferromagnetic exchange coupHng con- 
stant J = 1 in this work. The model has a SU{2) symme- 
try, which has been exploited in the DMRG program42ii^ 
We are interested in calculating the dynamical spin struc- 
ture factor 



can be expressed via 



S{k,uj) = /g (uj) 



- limIm(0|SJ;.- 



sjo) 



n rj-i-o ' ' '' UJ - H + b] + eo ^ 

TT ^—' 

n 

where |0) denotes the groundstate with energy eo- If we 
use open boundaries (obc) in the DMRG, we have to 
define the spin-operator Sfc in Fourier space (in the spirit 
of Ref . nil) as 



Sfc = 



-^'^^'^^WSj 



(2) 



j=i 



with quasi-momentum k = Itt/{L+ 1), I £N mod L. In 
this paper we also present results for periodic boundary 
conditions (pbc), even though the DMRG is not that 
efficient for pbc. The spin-operator is then given by 



Y^ exp{ijk)Sj 



(3) 



j=i 



with momentum k = It:/L, I G Z mod L. The asymp- 
totics of the dynamical structure factor are known from 
the analytical solutions,—"— and can be summarized as 
follows: 

• the dominant part of the structure factor comes 
from the two-spinon contributions. These are lo- 
cated in an interval bounded by 



cji = — I sinfcj and lu2 ~ 7r| sin — | 



(4) 



• the dynamical structure factor is known^>^ to di- 
verge as: 

S{k,Lu) oc (w-a;i)~Vln(l/(w-wi)) (5) 

for fc ^ TT, respectively 

5(7r,w) ccuj-'^.yinil/uj) (6) 

for k ~ n. 

III. CURRENT LANCZOS ALGORITHMS 

A. The basic Lanczos algorithm 

We are interested in calculating dynamic correlation 
functions at T = 0. In the Lehmann representation they 



Ia{uj) = -- limIm(0|A^- 



n i;^-o ' ' cj — H + 11] + eo 
= -y^\{E^\A\0)\^S{u-eo-e„). 



A |0) (7) 



Here \Ek) denote the eigenfunctions of the Hamiltonian 
H with eigenenergy e^. (spectral poles) . A is the operator 
that defines the correlation function. For small systems 
one can obtain the full set of eigenvectors and eigenen- 
ergies by exact diagonalization. However, this approach 
is limited rather quickly by the exponential growth of 
the Hilbert space. In order to circumvent this problem 
one can use an iterative eigensolver like for example the 
Lanczos algorithm. Within the Lanczos algorithm one 
recursively generates an orthogonal basis that spans the 
Krylov space. The recursion formula— for the so-called 
Lanczos vectors is given by 



\fn+l) = H\fn) - aMn) - bl\U-i) 
an^{fn\H\fn)/{fn\fn), 
bl^{fn\fn)/{.fn^l\.fn-l), ^0 = . 



(8) 



The full Hamiltonian is then mapped onto this Krylov 
space resulting in a tridiagonal effective Hamiltonian. 



H, 



off - 



/ao 6i 
6i a2 62 

62 ■■■ 



\ 







Vo 



'■■ '■• bn-1 
bn-1 an-ij 



(9) 



Because of the tridiagonal form this Hamiltonian can 
be easily diagonalized and it can be shown that the ex- 
tremal eigenvalues are good approximants^ to the ex- 
tremal eigenvalues of the original system. 

This procedure, originally developed for finding the ex- 
treme eigenvalues of large sparse Hamiltonians matrices, 
can now be adapted to the calculation of dynamical cor- 
relation functions^ For such functions as in (Eq. ([71)), 
we arc actually not interested in the full set of eigenvec- 
tors but only those vectors with non-zero overlap with 
the vector A\0). Therefore one starts the Lanczos algo- 
rithm with the vector |/o) = A|0)/(0|A^A|0), so that only 
eigenvectors and energies are created that give a non- 
zero contribution to the dynamical correlation function. 
If one diagonalizcs H^ff obtained in this way, its eigen- 
vectors and eigenvalues give direct access to the (lowest) 
poles and spectral weights of Eq. ([7]) . In order to obtain 
a "smooth" spectral function one can then broaden these 
delta functions. Alternatively, one can also calculate di- 
rectly a 'smooth' spectral function by using the continued 
fraction expansion 



/a(^) 



— lim Im 



(0|^t^|0) 



(10) 



z — ao 



with z = eo+uj + i-q. Both methods to obtain a 'smooth' 
spectral function are equivalent if one uses a Lorentz 
broadening of the delta peaks. The Lanczos method 
was used for our test case, the antiferromagnetic spin- 
^ Heisenberg chain (with periodic boundary conditions) 
by Fledderjohann et al^ and later by Karbach et al^ 
to obtain spectral weights and poles for chains up to 28 
sites. 

How well does this work? If applied to exact Hamil- 
tonians and states, the method is obviously limited by 
system size, as all "exact diagonalization" methods are. 
But there are more subtle issues. 

(i) The Lanczos algorithm will give the best conver- 
gence to the extremal eigenvectors that are contained in 
the starting vector. The interior eigenvectors will only 
converge after a large number of Lanczos iterations. This 
is in most cases not possible as numerical errors will de- 
stroy the orthogonality of the Lanczos vectors. This is the 
well-known "ghost problem"— i^ which leads to spurious 
(double) eigenvalues. Reorthogonalization and restarting 
methods are usually applied to tackle this problem. 

(ii) We are mainly interested in dynamical correlation 
functions in a continuous form (in the thermodynamic 
limit). In order to get a continuous form from our finite 
size discrete dynamical correlation functions we can use 
the two above described methods: Continued fraction ex- 
pansion or Lorentzian broadening of the spectral poles. 
For large systems and large broadenings this should give a 
decent approximation of the broadened correlation func- 
tion in the thermodynamic limit. Extracting the 77 — >■ 
limit a posteriori by removing the convolution necessar- 
ily requires at first a finite size scaling L — >■ 00 of the 
broadened data. Without a finite size scaling one would 
essentially try to 'deconvolute back' to a finite number 
of delta peaks. Furthermore the deconvolution is well- 
known to be a numerically ill-defined problem - there 
is no unique solution and the best solution can only be 
given by probability arguments (See Ref. |2J for details). 
In summary the extrapolation to the interesting thermo- 
dynamic limit 77 — > 0, i — > cx) is far from trivial especially 
as the two limits are not to be interchanged. 

This discussion is not specific to the Lanczos method, 
identical small shifts off the real frequency axis are 
mandatory in the correction vector method, too. In fact, 
the 77-broadening is also a feature of spectral functions 
obtained from time-dependent DMRG where the finite 
time-range of simulations is damped out by an artificial 
damping factor e~''* which leads to Lorentzian broaden- 
ing in frequency space. In the approach using Chebychev 
polynomials, damping, while taking a slightly different 
form, is also inevitable4£ 

(iii) Knowing the explicit positions of the spectral poles 
for a finite system is advantageous, because it allows, 
among other things, to perform proper finite size scal- 
ing for low-energy excitations and clearly identify possi- 
ble gaps in the spectrum. Furthermore, Holzner et al^ 
suggested recently to calculate the spectral weights for 
chains of different length and to use a systematic rescal- 



ing of the weights to approximate the dynamical corre- 
lation function in the thermodynamic limit. The idea 
of this rescaling is to simply translate the definition of 
the spectral function reading "spectral weight per unit 
frequency interval" into a mathematical formula4S We 
adapt their idea and use it for the spectral weights ex- 
tracted with the Lanczos method. We rescale each of the 
weights Vtn — |(-E„|Sfc|0)p by the width of the energy 
interval 



A„ 



{e^ 



n+l 



e„-i)/2 



(11) 



the pole is located in. For the lower bound cji of the 
spectral function this definition has to be replaced by 
Ai = (e2 -I- ei)/2 — wi. Using this rescaling scheme 
all spectral weights will now lie on the spectral func- 
tion calculated in the thermodynamic limit. Combining 
the spectral weights/poles of chains with different length, 
one can then generate a high resolution approximation to 
the spectral function in the thermodynamic limit and by 
this get a controlled approach to the problematic limit 
77 — > 0, L — >■ cxD. 



B. Lanczos with standard/adaptive DMRG 

Let us briefiy discuss how available Lanczos methods 
in the DMRG framework operate, regarding the broaden- 
ing and loss of orthogonality issues as well as additional, 
DMRG specific, approximations. 

The idea to use the Lanczos algorithm within standard 
DMRG was first put forward by Hallberg^ following the 
procedure outlined above, with |0) provided by a ground 
state DMRG calculation; operations H\'4}) are provided 
in DMRG such that the Lanczos algorithm can be im- 
plemented easily. Even if one assumes that |0) is avail- 
able at precision close to exact diagonalization, which 
is true to a very good approximation, an additional ap- 
proximation of this approach is that in the calculation 
of the Lanczos vectors (essentially i/"|/o)) H is avail- 
able only approximately and usually optimized for the 
operation H\Q). DMRG, if applied naively, will therefore 
introduce increasingly severe and systematic errors for 
the higher Lanczos vectors. A partial remedy is provided 
by the multi-targeting procedure, by which one finds a 
good approximate representation of several Lanczos vec- 
tors (which changes also the approximations made in the 
representation of H). The quandary one finds itself in 
is that if one targets few Lanczos vectors, those are rep- 
resented quite accurately, but all others quite badly; if 
one targets many Lanczos vectors, all are represented at 
similar, but quite low accuracy. Loss of orthogonality oc- 
curs as in every Lanczos procedure, but can be mended 
by reorthogonalization of the Lanczos vectors; as DMRG 
keeps them in a common basis representation, this can be 
done easily, but the price is the low accuracy of each indi- 
vidual Lanczos vector due to multiple targeting. Mostly 
the spectral function is presented as a smooth curve with 
an artificial broadening, but not as the set of poles and 



weights naturally following from the Lanczos representa- 
tion. This makes sense as the accuracy of individual pole 
positions and spectral weights is not so high in view of 
the systematic errors. 

Recently, some of the authors have proposed an adap- 
tive Lanczos methodic within a DMRG framework that is 
based on the same approach as Ref. |6|, but calculates the 
Lanczos vectors adaptively (i.e., using a multi-targeting 
approach, but restricting the Lanczos vectors to be tar- 
getted to the last three ones which occur in the recur- 
sion; this implies a continuous change of basis as the 
algorithm evolves). They could show that this approach 
is more efficient than the original formulation and allows 
one to directly calculate the position of the poles respec- 
tively the spectral weights for Green's function without 
any prior broadening because of higher accuracy. It must 
be stated however, that this approach is potentially vul- 
nerable to the ever-present problem of Lanczos "ghosts" , 
the appearance of spurious eigenvalues due to loss of or- 
thonormality between Lanczos vectors. The adaptive ba- 
sis changes make reorthogonalization schemes very diffi- 
cult if not effectively impossible. 



IV. LANCZOS WITH MATRIX PRODUCT 
STATES (MPS) 

Choosing the MPS formulation of DMRG in the Lanc- 
zos method will completely eliminate the need for mul- 
tiple targeting, giving additional accuracy, introduce an 
exact representation of H, and regain the possibility of a 
reorthogonalization scheme, lost in the adaptive Lanczos 
method. 

The algorithm proceeds exactly as before, the ground 
state |0) is written in MPS representation (as can be ob- 
tained from a standard ground state DMRG with minor 
modification or in a variational MPS approach). 

IV') = Yu A[ai]A[a2\..A[an-i\A[a^]\ai..(Jn) 



H is written as a matr ix p roduct operator (MPO); for 
examples, see Ref. [^ and [43. If we consider the recursion 
relation, we have to keep in mind that both the appli- 
cation of if to a state and the addition of two MPS (as 
occurs in the Lanczos procedure) lead, if carried out ex- 
actly, to a new MPS, but with a dimension that is larger 
than that of the original MPS: under addition, dimen- 
sions add; under application of an MPO to an MPS, the 
dimensions multiply (MPO dimensions are often 4 to 6). 
This is to be seen in contrast to the same operations in 
the conventional DMRG framework, where these opera- 
tions do not increase the numbers of DMRG states to be 
kept (i.e. the MPS matrix dimension), at the price of an 
approximate representation of the Hamiltonian H which 
is adapted only to specific Lanczos states and introduces 
quite severe inaccuracies for Lanczos states where H has 
been applied often. It is at this point that MPS based 



Lanczos dynamics has to introduce its only approxima- 
tion, namely the compression of the large-dimension MPS 
back to the original dimension at some loss of accuracy. 
This compression is carried out iteratively until the com- 
pressed state has become stable (for standard procedure, 
see Ref. |5|). This approximation avoids the systematic 
errors due to an approximate Hamiltonian and is in prac- 
tice less severe. Multiple targeting is eliminated. 

Within the MPS formulation we have found the best 
way to calculate the Lanczos states as MPS, but still the 
approximation due to the compression of the Lanczos 
states will lead to the ghost problem because the com- 
pression will increase the amount of orthogonality loss. In 
the MPS approach, this problem is not as bad as in the 
adaptive DMRG approachjiii but it is still exacerbated 
in comparison to the standard DMRG Lanczos method 
by Hallberg. The advantage is that this problem can be 
cured in comparison to the inexact Hamiltonian repre- 
sentation. There are many approaches to resolve the loss 
of orthogonality. The two most popular ones are total or 
partial re-orthogonalization and restarting proceduresi^ 
In a partial or total re-orthogonalization procedure one 
tries to re-orthogonalize the current Lanczos state to the 
previous ones by 

n — l _. 

l^n) = l/n)-E</"l^''')l'^')' 1^") =]^l^") (12) 
i 

after several steps or after each step. Here, A^ = 
((V'nlV'n))^ accounts for the proper normalization. We 
have tried partial and total re-orthogonalization (not 
shown), but it turns out that, due to the rather large 
overlap (/fc|/„) of the different Lanczos vectors, the cur- 
rent Lanczos vector will be changed significantly by this 
procedure; which results in an even worse representa- 
tion of the Hamiltonian. Alternatively speaking, the 
fact that re-orthogonalization involves sums implies again 
compression, which will make minor new approximations 
only if the re-orthogonalization is mild. Therefore we 
tried to devise a re-orthogonalization scheme which does 
not change the Lanczos vectors explicitly. 

To this end we rewrite the re-orthogonalization 
Eq. dni in the form 



IV'n) =YSin\fi) 



(13) 



where we have introduced the re-orthogonalization ma- 
trix Sin- It can be shown (see Appendix [X]) that Sin can 
be calculated from the overlap elements^ Wij = {fi\,fj) 
by a recursion relation. The matrix elements of the ef- 
fective Hamiltonian then become 



{lpn\H\iprn) 



n,7n 

E 



S*nSjm.{fi\H\fj) 



(14) 



By construction the vectors lipi) define an orthonormal 
basis system, i.e., on the left hand side we have ensured 
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FIG. 1. (Color online) Comparison of the convergence of the 
excitation energies and spectral weights for pseudomomen- 
tum fc = TT as function of Lanzcos iterations for the spin-i 
Heisenberg chain. The upper panel shows the behavior for 
the scheme without any reorthogonalization, the lower panel 
the new scheme. Calculations were done for a chain of length 
Z/ = 32 (obc), maximal MPS matrix dimension M = 512. 



that the effective Hamiltonian is connected by a unitary 
transformation to the right hand side. In particular, 
we do not have to change the calculated Lanczos vec- 
tors any more. The price that we are paying is that we 
need all overlap matrix elements Wij and matrix elements 
{fi\H\fj)r^ i.e., we have to store all Lanczos vectors 
generated during the calculation. Finally, the effective 
Hamiltonian given by the matrix elements {^lJn\H\■^pm) is 
not tridiagonal anymore and we have to apply a full di- 
agonalization to get the approximate excitation energies 
(The dimension of the effective Hamiltonian is given by 
the number of calculated Lanczos vectors and therefore 
the complete diagonalization is not computationally ex- 
pensive at all;so this additional cost is minor). 

Since the vectors {tpi) are orthonormal, the spectral 
weight flk belonging to the excitation energy e^. can now 



be calculated very easily fromSi 

Ofe/(0|iti|0) = \{Ek\A\0)\^ - |^Cfc,(V',;|Vo)l 



CfcO 



where we have used that the exact diagonalization of the 
effective Hamiltonian yields the (approximate) eigenvec- 
tors \Ek) = '^iCki\ipi)- In Fig. [1] a comparison of the 
convergence of the spectral weights of the MPS Lanczos 
method with (lower panel) and without (upper panel) 
re-orthogonalization is shown. One can clearly see that 
with a maximal MPS matrix dimension M = 512 only 
the first three spectral poles converge without the re- 
orthogonalization. With our reorthogonalization scheme 
one obtains a much better and more stable convergence, 
which allows to extract already the first 8 — 9 poles even 
for this small maximal MPS matrix dimension. Note 
however, that one of course does not get rid of all the 
ghosts. 

The obvious question within our algorithm thus is 
when to stop the Lanczos iterations and how to distin- 
guish real excitation energies from ghosts. The usual cri- 
terion for the quality of an (approximate) pole position 
uik is the residual. The residual vector is defined via 



rk 



H\Ek) - 0Jk\Ek) 



and the residual by 

Tk ^ {rkVk) = {Ek\{H - Ukf\Ek) . 

The residual can be viewed as a measure, how well the 
pole LOk and (normalized) vector \Ek) approximate an 
eigenvalue and eigenvector of H . We can again express 
this residual without knowing the eigenvectors explicitly, 
for the additional price that we have to measure all the 
matrix elements for H^ . Then the quantities entering the 
residual can be calculated as 

{E,\H\E,) = Y.cuCk, Y. S,nS,j{U\H\f,) , 

ij p,q 

{Ek\H^\Ek) = Jl^fe'^'^J^^'l^'l^j) 



— 2_^'^ki'^l'j 2_^ ^pi^qj\Jp\H \jq) 



With the help of the residuals we can define the follow- 
ing algorithm to calculate the position and weight of the 
poles: 

1. Remove all poles with a spectral weight below Vtcut 
and a residual that is larger than r^ut i because they 
very likely are ghosts. 

2. Follow the convergence of each of the remaining 
poles and take those with the smallest residual. 

In Fig. [21 the relative error in the position of the first 
three excitations for k = tt and the groundstate energy 
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FIG. 2. (Color online) The relative error of the first three 
spectral poles {k = 7r,pbc) calculated with the MPS Lanczos 
method in comparison to the two-spinon excitations from the 
Bethe Ansatz for different system sizes (upper panel) and for 
L = 32 for different maximal MPS matrix dimension M (lower 
panel). The groundstate (GS) error is also shown. 



calculated for periodic boundary conditions is shown. We 
compare our numerical results to exact ones based on the 
Bethe ansatz. The weight cutoff in our algorithm was 
chosen as Vtcut = 10~^ in all cases. The residual cutoff 
had to be increased with system size from r^ut = 0.1 for 
L = 24 to r^ut = 0.5 for L = 72 to be able to reliably 
extract the spectral weights. As expected one finds that 
the error in the groundstate energy is several orders of 
magnitude smaller than the error for the excitations. In 
the upper panel of Fig. [2] the system size L is varied 
for a constant maximal MPS matrix dimension M, while 
in the lower panel we keep L = 32 fixed and increase 
the maximal MPS matrix dimension M. The monotonic 
decrease of the error shows that it originates chiefly in 
the approximation of the Lanczos states by MPS. 

In the MPS-based Lanczos method we have now ob- 
tained accurate poles and spectral weights for a given 
finite chain length. In order to extrapolate these to the 
thermodynamic limit we will use the rcscaling scheme 



described in section UlI Al 



V. RESULTS 

1. k = TV 

We have tested our Lanczos algorithm for S{k = tt, uj) 
for both open and periodic boundary conditions. The 
results are collected in Fig. [31 The energy cutoff for 
the poles was chosen as i^cut = 10~^, and the residual 
cutoff had to be increased with system size to maximal 
fcut = 0.5 to be able to extract the spectral weights. 
One can clearly see that the discrete spectral weights, 
rescaled with the scheme discussed in section HIlAI nicely 
collapse onto a single curve, which is in good agreement 
with the results from the 2-spinon contributions of the 
Bethe ansatz in the thermodynamic limit,— independent 
of the chosen boundary conditions. Finite size effects 
are not very pronounced for fc = tt, and all (significant) 
spectral weight lies within the two-spinon bounds (see 
Eq. ^) 0Ji{k = tt) =0 to u}2{k = n) = n. The position 
of the second spectral pole differs substantially between 
open and periodic boundary conditions. Therefore for pe- 
riodic boundary conditions we would need much longer 
chains to fill the gap between the first and second spectral 
poles. The position of the first spectral pole moves very 
slowly with system size to the origin. Therefore, in order 
to resolve the divergence for w ^ we would have to go 
to much larger system sizes. The eigenstates that lead 
to non-vanishing weights for the dynamical spin struc- 
ture factor belong all to the S = 1 quantum sector i^ For 
fc = TT it is the eigenvector with the lowest energy in that 
subspace that gives the first spectral weight/pole. There- 
fore an alternative way to check the behavior for w — >■ 
for fc = TT is to do standard groundstate DMRG calcula- 
tions in that subspace (not shown). But in order to use 
our rcscaling scheme we need also the subsequent poles 
which are not as easy to obtain as long as translation 
symmetry is not implemented.— 

It turns out that the Lanczos method does not cap- 
ture the four-spinon weights (and higher spinon weights). 
In principle the Lanczos algorithm will also converge to- 
wards any multi-spinon state. But for this model the 
four-spinon and higher-spinon states show a rather strong 
finite size scaling,— which means that four-spinon states 
will appear in the low energy sector only for much longer 
chains than those considered here. Furthermore, these 
states will have a small weight in comparison to the two- 
spinon states and therefore they will be hard to distin- 
guish from ghost states or numerical noise. 

As mentioned earlier, the structure factor S{k = 7r,aj) 
shows a logarithmic correction to the leading divergence 
as w — >• (see Eq. ([S])). In order to resolve such an addi- 
tional feature it is instructive to plot ujS{k,uj,k), which is 
done in Fig. [5] for the periodic boundary conditions. Ap- 
parently, the data nicely fall onto the line from the Bethe 
ansatz, i.e. are in agreement with the predicted logarith- 
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FIG. 3. (Color online) Dynamical correlation function for 
k = n for open (upper panel) and periodic boundary (lower 
panel) conditions. The spectral weights/poles are chosen by 
the minimal residual scheme. 



FIG. 4. (Color online) Rescaled dynamical structure factor 
LjS{k,tjj) for k — TV and periodic boundary conditions. This 
figure is based on the same data as in Fig. |3l 
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mic divergence of this quantity. In order to make con- 
tact to earlier Lanczos methods, in Fig. [5l LuS{k = tt, lo) 
is plotted again, but for open boundary conditions and 
with an additional Lorentzian broadening of r; — 0.1. For 
this figure we do not have to use any residual cutoffs or 
weight cutoffs as for the extraction of the single weights. 
The broadening completely smears out the logarithmic 
divergence, but this is not specific to our method. In 
Figs. 8 and 14 of Ref. [la the same quantity up to a pro- 
portionality constant is discussed by Kiihner and White 
(their Fig. 8 is reproduced here for comparison. Fig. [6]). 
In their calculation, they used the standard DMRG Lanc- 
zos (broadened) and the (numerically much more ex- 
pensive) correction vector method. Although different 
schemes for dealing with the open boundary conditions 
preclude a quantitative comparison, their results clearly 
show that the standard Lanczos DMRG is not capable 
of producing a smooth, converged curve consistent with 
the predictions from Bethe ansatz, but shows substantial 
artificial structure. This comparison clearly shows that 



FIG. 5. (Color online) Rescaled dynamical structure factor 
ijjS{k,ij) for k = TT and open boundary conditions with a 
Lorentzian broadening of 77 = 0.1. 



wc have reached a definite improvement over the stan- 
dard Lanczos method. We want to emphasize that the 
inherent broadening of the correction vector completely 
annihilates any si gna ture of the divergence for a; — > 
(c.f. Fig. 14, Ref. Ha) . Note, however, that in order to 
unanibiguoulsy resolve or even predict such a logarith- 
mic divergence we would need to study system sizes well 
beyond anything possible presently. 



2. k = ^ 



In Fig.[7|we show S{k = f ,0;) (periodic boundary con- 
dictions) calculated with our method compared to the 2- 
spinon contribution obtained from Bethe ansatz 1^ The 
spectral signature starts now in the middle of the spec- 




FIG. 6. Dynamical structure factor ojS^{Lo,k) for k = n 
and open boundary conditions with a Lorentzian broadening 
of r; = 0.1, obtained in a standard Lanczos approach (with 
various numbers of target states Nl) and in the correction 
vector approach; reproduced from Fig. 8 of Ref. [la . Note 
that due to different schemes of dealing with the open bound- 
ary conditions results relate to those in Fig. [5] only up to a 
proportionality constant. 
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FIG. 7. (Color online) Dynamical correlation function for 
k = 7r/2 for periodic boundary conditions and M — 512. 
The spectral weights/poles are chosen by the minimal residual 
scheme. 



trum and therefore the extraction of poles and weights 
is much harder. We were able to obtain reasonable data 
up to a length L = 44 and M ~ 512. The low- lying exci- 
tations do not collapse as nicely onto the 2-spinon curve 
as for k = TT. One reason may lie in the problematic def- 
inition of the first interval in the rescaling scheme, which 
explicitly depends on the value of the lower bound, wi. 
It is evident that this value will also be subject to more 
or less strong finite size effects. These finite size effects 
are obviously even worse for open boundary conditions 



even though we can go to larger systems. 



VI. DISCUSSION 

As we have seen, the adaptation of the Lanczos 
algorithm for dynamical correlation functions to the 
framework of matrix product states produces much 
more accurate results than previous implementations 
in the DMRG framework and is capable of handling 
the broadly distributed spectral weight of the S = ^ 
Heisenberg antiferromagnet very well, which had been 
a major stumbling block for the original DMRG-based 
formulation. In more technical terms, the new algorithm 
offers three main advantages over previous Lanczos 
implementations in the DMRG: 

(i) We have shown that the Lanczos algorithm in 
MPS formulation is capable of calculating a number 
of discrete spectral weights/poles for spin chains that 
are longer than those that can be evaluated with exact 
diagonalization and the standard Lanczos algorithm. 
The algorithm can be easily implemented in any MPS 
algorithm as it just uses simple MPS algebra routines. It 
improves over the Lanczos implementations in classical 
DMRG and the adaptive Lanczos DMRG algorithm by 
avoiding multi-targeting and the inexact Hamiltonian 
representation. 

(ii) Furthermore the possibility of handling the Lanczos 
states individually by a single MPS allows for re- 
orthogonalization procedures that are not possible with 
the adaptive Lanczos method. This re-orthogonalization 
significantly improves the convergence of the spectral 
poles and makes the extraction of single weights much 
easier. 

(iii) As the number of discrete spectral poles is very 
limited (in the spin-^ Heisenberg chain) it is nearly 
impossible to get valuable information for the behavior 
in the thermodynamic limit from one single chain in 
this length scale. One should emphasize that all the 
information about the spectral function is just given by 
these discrete spectral weights/poles and therefore any 
method for the DMRG/MPS that gives a continuous 
function (e.g. correction vector, expansion in Chebyshev 
polynomials, Lanczos plus continuous fraction) just 
interpolates between these points depending on the 
broadening scheme of the used method (See section 
nil Ap . The used rescaling scheme based on the ex- 
plicit extraction of discrete spectral weights and poles 
gives a very controlled access to the problematic limit 
77 — > 0,ij — > cxD and by this gives more transparent 
results without broadening. 

We expect that the main competitor of this approach 
will be the other new numerical low-cost method 
provided by expansion in Chebychev polynomials. 
Compared to each other, both have advantages and 
disadvantages, such that no conclusive statement can 
be made at the moment, also in view of the lack of 



applications so far. Both rely on operations H\ip) 
carried out a similar number of times, but the Cheby- 
chev approach can handle longer chains (oc 300 sites) 
because of a surprising low maximal MPS matrix 
dimension (M ^ 32 — 64); this is arguably due to 
an inherent entanglement-reducing energy-truncation 
scheme. However, comparing directly the evaluation of 
discrete spectral weights/poles, one has to perform many 
iterations (500-1000) within the Chebyshev expansion. 
Furthermore, one has to fit the broadened spectral 
weights to get the discrete poles/weights including a fit 
error. This becomes even more challenging for dense 
spectral poles, where the broadened peaks overlap. Gen- 
erally, broadening is an inherent part of the Chebychev 
approach. It also has more tuning parameters that have 
to be optimized for the method to work properly. 

Some algorithms for time evolution of quantum states 
also use a Krylov space representation in Matrix Product 
States 1^1^"— Here our proposed re-orthogonalization 
may also be useful. 

Recently, two algorithms for calculating spectra of 
translationally invariant Hamiltonians with (unitary) 
MPS working directly in the thermodynamic limit^/ 
with periodic boundary conditions^ were suggested 
that can approximate many branches of the dispersion 
relation with high precision. It would be really interest- 
ing to see whether these algorithms can also calculate 
dynamic correlation functions with that precision and 
how they compare to the Lanczos expansion. 
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Appendix A: Recursive formula for 
reorthogonalization 

We want to express the reorthogonalization equation 

n— 1 -. 

\^n)^\.fn)-Y.^In\i^^)m, IV'n) = ^Hn) (Al) 

by a matr ix S with j-iAn) = YJi Sin\Ii) (with iV„ = 
(V'nlV'n))- This leads to: 

n-l 

■i 

n — 1 i 

= I/") - ^^hnSk'i\fk'), 
i k' 
i 
kin = E SkiWnk, Wij = {f^\fj) . 



Now one can deduce a recursion formula for the matrix 
elements of S: 



-'pn 



Z^q=p kqnbpq, 11 p < Tl 

1, a p = n 



(A2) 



with 



5, 



pn- ^^i>pn, ^00-^^ 



Nn 



. E SpnSqnWpq, No=VWi 

\ g,p=o 



00 • 
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